Apply WGCNA to Inv4m DEGs to identify co-expression modules.
Approach: Unsigned networks with power = 12, minModuleSize = 10, better suited for identifying regulatory hubs like JMJ genes that can both activate and repress targets.
library(SummarizedExperiment)
library(tidyverse)
library(WGCNA)
library(pheatmap)
voomR <- readRDS(file.path(paths$intermediate, "normalized_expression_voom_object_leaf_trt.rds"))
cat("Expression data loaded:\n")
## Expression data loaded:
cat(" Genes:", nrow(voomR$E), "\n")
## Genes: 24249
cat(" Samples:", ncol(voomR$E), "\n")
## Samples: 43
effects_df <- read_csv(
file.path(paths$intermediate, "predictor_effects_leaf_interaction_model.csv")
)
inv4m_degs <- effects_df %>%
filter(
predictor == "Inv4m",
is_DEG == TRUE
) %>%
pull(gene) %>%
unique()
cat("Inv4m DEGs identified:", length(inv4m_degs), "\n")
## Inv4m DEGs identified: 776
cat(" Upregulated:",
sum(effects_df$predictor == "Inv4m" &
effects_df$is_DEG &
effects_df$regulation == "Upregulated"), "\n")
## Upregulated: 237
cat(" Downregulated:",
sum(effects_df$predictor == "Inv4m" &
effects_df$is_DEG &
effects_df$regulation == "Downregulated"), "\n")
## Downregulated: 539
expr_matrix_degs <- t(voomR$E[inv4m_degs, ])
cat("Expression matrix dimensions:", dim(expr_matrix_degs), "\n")
## Expression matrix dimensions: 43 776
cat(" Samples (rows):", nrow(expr_matrix_degs), "\n")
## Samples (rows): 43
cat(" Genes (columns):", ncol(expr_matrix_degs), "\n")
## Genes (columns): 776
sample_table <- voomR$targets %>%
as.data.frame() %>%
rownames_to_column("Sample") %>%
select(
Sample,
Genotype,
leaf_tissue,
Treatment
) %>%
mutate(
Genotype = factor(Genotype, levels = c("CTRL", "Inv4m")),
leaf_tissue = factor(leaf_tissue)
)
rownames(sample_table) <- sample_table$Sample
cat("Sample distribution:\n")
## Sample distribution:
with(sample_table, {
cat("Genotype × Leaf tissue:\n")
print(table(Genotype, leaf_tissue))
cat("\nGenotype × Treatment:\n")
print(table(Genotype, Treatment))
})
## Genotype × Leaf tissue:
## leaf_tissue
## Genotype 1 2 3 4
## CTRL 4 6 5 5
## Inv4m 7 5 6 5
##
## Genotype × Treatment:
## Treatment
## Genotype +P -P
## CTRL 12 8
## Inv4m 12 11
datExpr <- expr_matrix_degs
cat("Expression data prepared for WGCNA:\n")
## Expression data prepared for WGCNA:
cat(" Samples (rows):", nrow(datExpr), "\n")
## Samples (rows): 43
cat(" Genes (columns):", ncol(datExpr), "\n")
## Genes (columns): 776
powers <- c(1:20)
sft_all <- pickSoftThreshold(
datExpr,
powerVector = powers,
networkType = "unsigned",
verbose = 5
)
## pickSoftThreshold: will use block size 776.
## pickSoftThreshold: calculating connectivity for given powers...
## ..working on genes 1 through 776 of 776
## Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
## 1 1 0.49400 2.6700 0.936 324.00 314.000 468.0
## 2 2 0.23300 0.8980 0.813 172.00 161.000 309.0
## 3 3 0.00295 -0.0599 0.672 102.00 93.000 222.0
## 4 4 0.21500 -0.4490 0.671 66.10 58.000 168.0
## 5 5 0.58600 -0.7990 0.716 45.20 37.500 134.0
## 6 6 0.76900 -1.0200 0.792 32.40 25.300 112.0
## 7 7 0.86000 -1.1400 0.855 24.10 17.600 96.4
## 8 8 0.91700 -1.2000 0.902 18.40 12.600 84.5
## 9 9 0.93200 -1.2400 0.915 14.50 8.920 75.2
## 10 10 0.93200 -1.2500 0.912 11.60 6.490 67.7
## 11 11 0.95700 -1.2700 0.946 9.53 4.810 61.5
## 12 12 0.93200 -1.2700 0.915 7.91 3.580 56.3
## 13 13 0.94400 -1.2600 0.935 6.66 2.690 51.8
## 14 14 0.94100 -1.2600 0.929 5.67 2.070 47.9
## 15 15 0.96600 -1.2400 0.961 4.88 1.560 44.5
## 16 16 0.96700 -1.2200 0.964 4.23 1.230 41.5
## 17 17 0.95600 -1.2300 0.950 3.70 0.936 38.7
## 18 18 0.96200 -1.2100 0.961 3.26 0.739 36.3
## 19 19 0.95700 -1.2100 0.951 2.89 0.581 34.1
## 20 20 0.96700 -1.1900 0.963 2.57 0.465 32.1
par(mfrow = c(1, 2))
plot(
sft_all$fitIndices[, 1],
-sign(sft_all$fitIndices[, 3]) * sft_all$fitIndices[, 2],
xlab = "Soft Threshold (power)",
ylab = "Scale Free Topology Model Fit, signed R^2",
main = "Scale Independence - All Samples",
type = "n"
)
text(
sft_all$fitIndices[, 1],
-sign(sft_all$fitIndices[, 3]) * sft_all$fitIndices[, 2],
labels = powers,
cex = 0.9,
col = "red"
)
abline(h = 0.80, col = "blue")
plot(
sft_all$fitIndices[, 1],
sft_all$fitIndices[, 5],
xlab = "Soft Threshold (power)",
ylab = "Mean Connectivity",
main = "Mean Connectivity - All Samples",
type = "n"
)
text(
sft_all$fitIndices[, 1],
sft_all$fitIndices[, 5],
labels = powers,
cex = 0.9,
col = "red"
)
pow_all <- sft_all$powerEstimate
if (is.na(pow_all)) {
pow_all <- sft_all$fitIndices %>%
filter(-sign(Slope) * SFT.R.sq > 0.80) %>%
slice_head(n = 1) %>%
pull(Power)
if (length(pow_all) == 0) {
pow_all <- 12
cat("\nWarning: No power achieved R^2 > 0.80, using default power = 12\n")
}
}
cat("\nSelected soft-thresholding power (all samples):", pow_all, "\n")
##
## Selected soft-thresholding power (all samples): 7
cat("Scale-free topology fit R^2:",
sft_all$fitIndices[sft_all$fitIndices$Power == pow_all, "SFT.R.sq"], "\n")
## Scale-free topology fit R^2: 0.8601684
cat("Building UNSIGNED network with blockwiseModules:\n")
## Building UNSIGNED network with blockwiseModules:
net_all <- blockwiseModules(
datExpr,
power = pow_all,
networkType = "unsigned",
TOMType = "unsigned",
minModuleSize = 10,
maxBlockSize = 25000,
reassignThreshold = 0,
mergeCutHeight = 0,
numericLabels = TRUE,
pamRespectsDendro = FALSE,
deepSplit = 4,
verbose = 3
)
## Calculating module eigengenes block-wise from all genes
## Flagging genes and samples with too many missing values...
## ..step 1
## ..Working on block 1 .
## TOM calculation: adjacency..
## ..will not use multithreading.
## Fraction of slow calculations: 0.000000
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.
## ....clustering..
## ....detecting modules..
## ....calculating module eigengenes..
## ....checking kME in modules..
## ..merging modules that are too close..
## mergeCloseModules: Merging modules whose distance is less than 0
## Calculating new MEs...
wgcna_all <- list(
combined = structure(
list(
datExpr = datExpr,
moduleColors = net_all$colors,
MEs = net_all$MEs,
sampleTable = sample_table
),
class = "WGCNA"
)
)
cat("\nNetwork built successfully\n")
##
## Network built successfully
cat(" Modules:", length(unique(net_all$colors)), "\n")
## Modules: 14
cat(" Genes:", length(net_all$colors), "\n")
## Genes: 776
ctrl_samples <- rownames(sample_table)[sample_table$Genotype == "CTRL"]
inv4m_samples <- rownames(sample_table)[sample_table$Genotype == "Inv4m"]
datExpr_ctrl <- datExpr[ctrl_samples, ]
datExpr_inv4m <- datExpr[inv4m_samples, ]
cat("CTRL expression data:\n")
## CTRL expression data:
cat(" Samples:", nrow(datExpr_ctrl), "\n")
## Samples: 20
cat(" Genes:", ncol(datExpr_ctrl), "\n")
## Genes: 776
cat("\nInv4m expression data:\n")
##
## Inv4m expression data:
cat(" Samples:", nrow(datExpr_inv4m), "\n")
## Samples: 23
cat(" Genes:", ncol(datExpr_inv4m), "\n")
## Genes: 776
sft_ctrl <- pickSoftThreshold(
datExpr_ctrl,
powerVector = powers,
networkType = "unsigned",
verbose = 5
)
## pickSoftThreshold: will use block size 776.
## pickSoftThreshold: calculating connectivity for given powers...
## ..working on genes 1 through 776 of 776
## Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
## 1 1 0.4840 1.140 0.944 254.00 256.000 373.0
## 2 2 0.0303 -0.125 0.832 121.00 116.000 235.0
## 3 3 0.4830 -0.555 0.849 68.40 60.900 166.0
## 4 4 0.7670 -0.764 0.930 42.90 35.500 124.0
## 5 5 0.8230 -0.874 0.927 28.80 21.700 96.4
## 6 6 0.8610 -0.949 0.929 20.40 13.800 76.9
## 7 7 0.9070 -0.985 0.953 14.90 9.020 62.5
## 8 8 0.9250 -1.020 0.962 11.30 6.320 51.6
## 9 9 0.9500 -1.060 0.968 8.72 4.430 43.3
## 10 10 0.9400 -1.070 0.960 6.87 3.170 37.0
## 11 11 0.9420 -1.110 0.962 5.50 2.340 32.0
## 12 12 0.9510 -1.120 0.976 4.47 1.720 27.9
## 13 13 0.9450 -1.150 0.971 3.67 1.310 24.5
## 14 14 0.9520 -1.160 0.982 3.05 1.010 21.7
## 15 15 0.9540 -1.180 0.979 2.56 0.781 19.2
## 16 16 0.9250 -1.210 0.956 2.16 0.612 17.1
## 17 17 0.9320 -1.220 0.964 1.84 0.474 15.3
## 18 18 0.9420 -1.220 0.977 1.58 0.369 13.8
## 19 19 0.9450 -1.230 0.977 1.37 0.294 12.4
## 20 20 0.9290 -1.240 0.948 1.19 0.232 11.2
sft_inv4m <- pickSoftThreshold(
datExpr_inv4m,
powerVector = powers,
networkType = "unsigned",
verbose = 5
)
## pickSoftThreshold: will use block size 776.
## pickSoftThreshold: calculating connectivity for given powers...
## ..working on genes 1 through 776 of 776
## Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
## 1 1 0.5030 2.060 0.881 225.000 230.0000 310.00
## 2 2 0.0618 0.329 0.749 95.400 97.1000 167.00
## 3 3 0.0678 -0.277 0.749 48.500 47.7000 103.00
## 4 4 0.3080 -0.681 0.782 27.500 25.9000 68.60
## 5 5 0.4650 -0.870 0.842 16.800 15.0000 48.20
## 6 6 0.5940 -1.000 0.896 10.900 9.1800 35.20
## 7 7 0.6710 -1.080 0.932 7.370 5.8300 26.60
## 8 8 0.7240 -1.170 0.938 5.160 3.8400 20.50
## 9 9 0.7650 -1.190 0.947 3.720 2.6000 16.20
## 10 10 0.7730 -1.250 0.931 2.740 1.7800 12.90
## 11 11 0.8030 -1.260 0.943 2.070 1.2500 10.50
## 12 12 0.8040 -1.300 0.925 1.590 0.8930 8.66
## 13 13 0.8290 -1.310 0.938 1.240 0.6550 7.21
## 14 14 0.8480 -1.310 0.935 0.981 0.4800 6.06
## 15 15 0.8860 -1.280 0.942 0.786 0.3660 5.14
## 16 16 0.8240 -1.330 0.855 0.638 0.2680 4.48
## 17 17 0.8920 -1.410 0.955 0.523 0.2000 4.34
## 18 18 0.2720 -2.640 0.208 0.434 0.1500 4.22
## 19 19 0.3050 -3.650 0.257 0.363 0.1180 4.12
## 20 20 0.3200 -3.670 0.275 0.306 0.0889 4.02
par(mfrow = c(2, 2))
plot(
sft_ctrl$fitIndices[, 1],
-sign(sft_ctrl$fitIndices[, 3]) * sft_ctrl$fitIndices[, 2],
xlab = "Soft Threshold (power)",
ylab = "Scale Free Topology Fit R^2",
main = "Scale Independence - CTRL",
type = "n"
)
text(
sft_ctrl$fitIndices[, 1],
-sign(sft_ctrl$fitIndices[, 3]) * sft_ctrl$fitIndices[, 2],
labels = powers,
cex = 0.9,
col = "red"
)
abline(h = 0.80, col = "blue")
plot(
sft_ctrl$fitIndices[, 1],
sft_ctrl$fitIndices[, 5],
xlab = "Soft Threshold (power)",
ylab = "Mean Connectivity",
main = "Mean Connectivity - CTRL",
type = "n"
)
text(
sft_ctrl$fitIndices[, 1],
sft_ctrl$fitIndices[, 5],
labels = powers,
cex = 0.9,
col = "red"
)
plot(
sft_inv4m$fitIndices[, 1],
-sign(sft_inv4m$fitIndices[, 3]) * sft_inv4m$fitIndices[, 2],
xlab = "Soft Threshold (power)",
ylab = "Scale Free Topology Fit R^2",
main = "Scale Independence - Inv4m",
type = "n"
)
text(
sft_inv4m$fitIndices[, 1],
-sign(sft_inv4m$fitIndices[, 3]) * sft_inv4m$fitIndices[, 2],
labels = powers,
cex = 0.9,
col = "red"
)
abline(h = 0.80, col = "blue")
plot(
sft_inv4m$fitIndices[, 1],
sft_inv4m$fitIndices[, 5],
xlab = "Soft Threshold (power)",
ylab = "Mean Connectivity",
main = "Mean Connectivity - Inv4m",
type = "n"
)
text(
sft_inv4m$fitIndices[, 1],
sft_inv4m$fitIndices[, 5],
labels = powers,
cex = 0.9,
col = "red"
)
pow_ctrl <- sft_ctrl$powerEstimate
if (is.na(pow_ctrl)) {
pow_ctrl <- sft_ctrl$fitIndices %>%
filter(-sign(Slope) * SFT.R.sq > 0.80) %>%
slice_head(n = 1) %>%
pull(Power)
if (length(pow_ctrl) == 0) {
pow_ctrl <- 12
cat("\nWarning: CTRL - No power achieved R^2 > 0.80, using default = 12\n")
}
}
pow_inv4m <- sft_inv4m$powerEstimate
if (is.na(pow_inv4m)) {
pow_inv4m <- sft_inv4m$fitIndices %>%
filter(-sign(Slope) * SFT.R.sq > 0.80) %>%
slice_head(n = 1) %>%
pull(Power)
if (length(pow_inv4m) == 0) {
pow_inv4m <- 12
cat("\nWarning: Inv4m - No power achieved R^2 > 0.80, using default = 12\n")
}
}
cat("\nSelected powers:\n")
##
## Selected powers:
cat(" CTRL:", pow_ctrl,
"(R^2 =", sft_ctrl$fitIndices[sft_ctrl$fitIndices$Power == pow_ctrl, "SFT.R.sq"], ")\n")
## CTRL: 6 (R^2 = 0.8605398 )
cat(" Inv4m:", pow_inv4m,
"(R^2 =", sft_inv4m$fitIndices[sft_inv4m$fitIndices$Power == pow_inv4m, "SFT.R.sq"], ")\n")
## Inv4m: 15 (R^2 = 0.8863985 )
cat("Building CTRL network:\n")
## Building CTRL network:
net_ctrl <- blockwiseModules(
datExpr_ctrl,
power = pow_ctrl,
networkType = "unsigned",
TOMType = "unsigned",
minModuleSize = 10,
maxBlockSize = 25000,
reassignThreshold = 0,
mergeCutHeight = 0,
numericLabels = TRUE,
pamRespectsDendro = FALSE,
deepSplit = 4,
verbose = 3
)
## Calculating module eigengenes block-wise from all genes
## Flagging genes and samples with too many missing values...
## ..step 1
## ..Working on block 1 .
## TOM calculation: adjacency..
## ..will not use multithreading.
## Fraction of slow calculations: 0.000000
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.
## ....clustering..
## ....detecting modules..
## ....calculating module eigengenes..
## ....checking kME in modules..
## ..removing 6 genes from module 1 because their KME is too low.
## ..removing 11 genes from module 2 because their KME is too low.
## ..removing 1 genes from module 3 because their KME is too low.
## ..removing 1 genes from module 4 because their KME is too low.
## ..removing 3 genes from module 8 because their KME is too low.
## ..removing 3 genes from module 10 because their KME is too low.
## ..removing 1 genes from module 15 because their KME is too low.
## ..merging modules that are too close..
## mergeCloseModules: Merging modules whose distance is less than 0
## Calculating new MEs...
wgcna_ctrl <- list(
ctrl = structure(
list(
datExpr = datExpr_ctrl,
moduleColors = net_ctrl$colors,
MEs = net_ctrl$MEs,
sampleTable = sample_table[ctrl_samples, ]
),
class = "WGCNA"
)
)
cat("\nCTRL network built\n")
##
## CTRL network built
cat(" Modules:", length(unique(net_ctrl$colors)), "\n")
## Modules: 19
cat(" Genes:", length(net_ctrl$colors), "\n")
## Genes: 776
cat("Building Inv4m network:\n")
## Building Inv4m network:
net_inv4m <- blockwiseModules(
datExpr_inv4m,
power = pow_inv4m,
networkType = "unsigned",
TOMType = "unsigned",
minModuleSize = 10,
maxBlockSize = 25000,
reassignThreshold = 0,
mergeCutHeight = 0,
numericLabels = TRUE,
pamRespectsDendro = FALSE,
deepSplit = 4,
verbose = 3
)
## Calculating module eigengenes block-wise from all genes
## Flagging genes and samples with too many missing values...
## ..step 1
## ..Working on block 1 .
## TOM calculation: adjacency..
## ..will not use multithreading.
## Fraction of slow calculations: 0.000000
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.
## ....clustering..
## ....detecting modules..
## ....calculating module eigengenes..
## ....checking kME in modules..
## ..merging modules that are too close..
## mergeCloseModules: Merging modules whose distance is less than 0
## Calculating new MEs...
wgcna_inv4m <- list(
inv4m = structure(
list(
datExpr = datExpr_inv4m,
moduleColors = net_inv4m$colors,
MEs = net_inv4m$MEs,
sampleTable = sample_table[inv4m_samples, ]
),
class = "WGCNA"
)
)
cat("\nInv4m network built\n")
##
## Inv4m network built
cat(" Modules:", length(unique(net_inv4m$colors)), "\n")
## Modules: 19
cat(" Genes:", length(net_inv4m$colors), "\n")
## Genes: 776
n_modules_all <- length(unique(net_all$colors))
module_sizes_all <- table(net_all$colors)
cat("\nCombined network (all samples):\n")
##
## Combined network (all samples):
cat(" Total modules:", n_modules_all, "\n")
## Total modules: 14
cat(" Total genes:", length(net_all$colors), "\n")
## Total genes: 776
cat("\nModule sizes:\n")
##
## Module sizes:
print(sort(module_sizes_all, decreasing = TRUE))
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14
## 340 77 60 50 35 31 31 30 29 29 21 20 13 10
saveRDS(wgcna_all, file.path(paths$intermediate, "wgcna_all.rds"))
cat("\nNetwork saved to wgcna_all.rds\n")
##
## Network saved to wgcna_all.rds
n_modules_ctrl <- length(unique(net_ctrl$colors))
module_sizes_ctrl <- table(net_ctrl$colors)
n_modules_inv4m <- length(unique(net_inv4m$colors))
module_sizes_inv4m <- table(net_inv4m$colors)
cat("\nCTRL network:\n")
##
## CTRL network:
cat(" Total modules:", n_modules_ctrl, "\n")
## Total modules: 19
cat(" Total genes:", length(net_ctrl$colors), "\n")
## Total genes: 776
cat("\nModule sizes:\n")
##
## Module sizes:
print(sort(module_sizes_ctrl, decreasing = TRUE))
##
## 1 2 3 4 5 6 7 8 9 10 11 0 12 13 14 15 16 17 18
## 116 102 82 63 56 40 39 36 36 30 28 27 27 18 17 16 15 15 13
cat("\n\nInv4m network:\n")
##
##
## Inv4m network:
cat(" Total modules:", n_modules_inv4m, "\n")
## Total modules: 19
cat(" Total genes:", length(net_inv4m$colors), "\n")
## Total genes: 776
cat("\nModule sizes:\n")
##
## Module sizes:
print(sort(module_sizes_inv4m, decreasing = TRUE))
##
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
## 288 62 50 45 44 40 38 33 26 21 19 17 16 15 14 13 12 12 11
saveRDS(wgcna_ctrl, file.path(paths$intermediate, "wgcna_ctrl.rds"))
saveRDS(wgcna_inv4m, file.path(paths$intermediate, "wgcna_inv4m.rds"))
cat("\nGenotype-specific networks saved\n")
##
## Genotype-specific networks saved
MEs_all <- net_all$MEs
ME_data_all <- MEs_all %>%
as.data.frame() %>%
rownames_to_column("Sample") %>%
left_join(
sample_table %>% select(Sample, Genotype, leaf_tissue, Treatment),
by = "Sample"
)
module_names_all <- grep("^ME", colnames(MEs_all), value = TRUE)
dme_results <- map_dfr(module_names_all, function(mod) {
me_values <- ME_data_all[[mod]]
genotype <- ME_data_all$Genotype
inv4m_vals <- me_values[genotype == "Inv4m"]
ctrl_vals <- me_values[genotype == "CTRL"]
t_result <- t.test(inv4m_vals, ctrl_vals)
tibble(
module = mod,
mean_Inv4m = mean(inv4m_vals),
mean_CTRL = mean(ctrl_vals),
diff = mean_Inv4m - mean_CTRL,
t_statistic = t_result$statistic,
p_value = t_result$p.value
)
})
dme_results <- dme_results %>%
mutate(fdr = p.adjust(p_value, method = "fdr")) %>%
arrange(p_value)
cat("Differential module expression results:\n")
## Differential module expression results:
print(dme_results)
## # A tibble: 14 × 7
## module mean_Inv4m mean_CTRL diff t_statistic p_value fdr
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ME1 -0.141 0.162 -0.303 -51.5 7.30e-39 1.02e-37
## 2 ME5 -0.100 0.115 -0.215 -6.33 1.61e- 7 1.13e- 6
## 3 ME9 -0.0988 0.114 -0.212 -6.08 5.76e- 7 2.69e- 6
## 4 ME7 -0.0909 0.104 -0.195 -5.28 5.39e- 6 1.70e- 5
## 5 ME6 0.0935 -0.107 0.201 5.42 6.08e- 6 1.70e- 5
## 6 ME13 -0.0838 0.0964 -0.180 -4.69 3.11e- 5 7.27e- 5
## 7 ME14 -0.0834 0.0959 -0.179 -4.60 4.32e- 5 8.65e- 5
## 8 ME10 -0.0821 0.0944 -0.177 -4.45 8.25e- 5 1.44e- 4
## 9 ME8 0.0778 -0.0894 0.167 4.24 1.25e- 4 1.95e- 4
## 10 ME3 -0.0774 0.0890 -0.166 -4.04 3.16e- 4 4.42e- 4
## 11 ME11 -0.0737 0.0848 -0.159 -3.77 6.82e- 4 8.68e- 4
## 12 ME4 -0.0606 0.0697 -0.130 -2.92 6.49e- 3 7.19e- 3
## 13 ME2 -0.0597 0.0686 -0.128 -2.89 6.68e- 3 7.19e- 3
## 14 ME12 0.0563 -0.0648 0.121 2.71 1.04e- 2 1.04e- 2
sig_modules <- dme_results %>%
filter(fdr < 0.05)
cat("\nSignificant genotype-associated modules (FDR < 0.05):",
nrow(sig_modules), "\n")
##
## Significant genotype-associated modules (FDR < 0.05): 14
if (nrow(sig_modules) > 0) {
print(sig_modules)
}
## # A tibble: 14 × 7
## module mean_Inv4m mean_CTRL diff t_statistic p_value fdr
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ME1 -0.141 0.162 -0.303 -51.5 7.30e-39 1.02e-37
## 2 ME5 -0.100 0.115 -0.215 -6.33 1.61e- 7 1.13e- 6
## 3 ME9 -0.0988 0.114 -0.212 -6.08 5.76e- 7 2.69e- 6
## 4 ME7 -0.0909 0.104 -0.195 -5.28 5.39e- 6 1.70e- 5
## 5 ME6 0.0935 -0.107 0.201 5.42 6.08e- 6 1.70e- 5
## 6 ME13 -0.0838 0.0964 -0.180 -4.69 3.11e- 5 7.27e- 5
## 7 ME14 -0.0834 0.0959 -0.179 -4.60 4.32e- 5 8.65e- 5
## 8 ME10 -0.0821 0.0944 -0.177 -4.45 8.25e- 5 1.44e- 4
## 9 ME8 0.0778 -0.0894 0.167 4.24 1.25e- 4 1.95e- 4
## 10 ME3 -0.0774 0.0890 -0.166 -4.04 3.16e- 4 4.42e- 4
## 11 ME11 -0.0737 0.0848 -0.159 -3.77 6.82e- 4 8.68e- 4
## 12 ME4 -0.0606 0.0697 -0.130 -2.92 6.49e- 3 7.19e- 3
## 13 ME2 -0.0597 0.0686 -0.128 -2.89 6.68e- 3 7.19e- 3
## 14 ME12 0.0563 -0.0648 0.121 2.71 1.04e- 2 1.04e- 2
ME_heatmap_data <- ME_data_all %>%
select(Sample, Genotype, leaf_tissue, Treatment, all_of(module_names_all))
anno_df <- ME_heatmap_data %>%
select(Sample, Genotype, leaf_tissue, Treatment) %>%
column_to_rownames("Sample")
ME_matrix <- ME_heatmap_data %>%
select(all_of(module_names_all)) %>%
as.matrix()
rownames(ME_matrix) <- ME_heatmap_data$Sample
pheatmap(
t(ME_matrix),
annotation_col = anno_df,
scale = "row",
clustering_distance_rows = "correlation",
clustering_distance_cols = "euclidean",
main = "Module Eigengenes by Sample",
fontsize = 10
)
sig_module_plots <- sig_modules %>%
pull(module) %>%
head(4)
if (length(sig_module_plots) > 0) {
plot_data <- ME_data_all %>%
select(Sample, Genotype, leaf_tissue, all_of(sig_module_plots)) %>%
pivot_longer(
cols = all_of(sig_module_plots),
names_to = "module",
values_to = "eigengene"
)
ggplot(plot_data, aes(x = Genotype, y = eigengene, fill = Genotype)) +
geom_boxplot() +
geom_jitter(width = 0.2, alpha = 0.5) +
facet_wrap(~ module, scales = "free_y") +
theme_classic(base_size = 12) +
labs(
title = "Top Genotype-Associated Modules",
y = "Module Eigengene",
x = NULL
) +
theme(legend.position = "none")
}
mod_all <- data.frame(
gene = names(net_all$colors),
module = paste0("ME", net_all$colors)
)
module_deg_stats_all <- mod_all %>%
left_join(
effects_df %>%
filter(predictor == "Inv4m") %>%
select(gene, logFC, adj.P.Val, regulation, in_Inv4m, in_cis),
by = "gene"
)
cat("Module assignments (all samples) with DEG statistics:\n")
## Module assignments (all samples) with DEG statistics:
cat(" Total genes:", nrow(module_deg_stats_all), "\n")
## Total genes: 776
cat(" Genes with stats:", sum(!is.na(module_deg_stats_all$logFC)), "\n")
## Genes with stats: 776
mod_ctrl <- data.frame(
gene = names(net_ctrl$colors),
module = paste0("ME", net_ctrl$colors)
)
module_deg_stats_ctrl <- mod_ctrl %>%
left_join(
effects_df %>%
filter(predictor == "Inv4m") %>%
select(gene, logFC, adj.P.Val, regulation, in_Inv4m, in_cis),
by = "gene"
)
mod_inv4m <- data.frame(
gene = names(net_inv4m$colors),
module = paste0("ME", net_inv4m$colors)
)
module_deg_stats_inv4m <- mod_inv4m %>%
left_join(
effects_df %>%
filter(predictor == "Inv4m") %>%
select(gene, logFC, adj.P.Val, regulation, in_Inv4m, in_cis),
by = "gene"
)
cat("\nModule assignments (CTRL) with DEG statistics:\n")
##
## Module assignments (CTRL) with DEG statistics:
cat(" Total genes:", nrow(module_deg_stats_ctrl), "\n")
## Total genes: 776
cat(" Genes with stats:", sum(!is.na(module_deg_stats_ctrl$logFC)), "\n")
## Genes with stats: 776
cat("\nModule assignments (Inv4m) with DEG statistics:\n")
##
## Module assignments (Inv4m) with DEG statistics:
cat(" Total genes:", nrow(module_deg_stats_inv4m), "\n")
## Total genes: 776
cat(" Genes with stats:", sum(!is.na(module_deg_stats_inv4m$logFC)), "\n")
## Genes with stats: 776
cat("\nCalculating adjacency matrix (all samples)...\n")
##
## Calculating adjacency matrix (all samples)...
adj_all <- adjacency(
datExpr,
power = pow_all,
type = "unsigned"
)
cat("Adjacency matrix dimensions:", dim(adj_all), "\n")
## Adjacency matrix dimensions: 776 776
edge_indices_all <- which(adj_all > 0 & upper.tri(adj_all, diag = FALSE),
arr.ind = TRUE)
edge_list_all <- data.frame(
from = rownames(adj_all)[edge_indices_all[, 1]],
to = colnames(adj_all)[edge_indices_all[, 2]],
weight = adj_all[edge_indices_all]
) %>%
arrange(desc(weight))
cat("Total edges (weight > 0):", nrow(edge_list_all), "\n")
## Total edges (weight > 0): 300700
cat("\nCalculating adjacency matrices (genotype-specific)...\n")
##
## Calculating adjacency matrices (genotype-specific)...
adj_ctrl <- adjacency(
datExpr_ctrl,
power = pow_ctrl,
type = "unsigned"
)
adj_inv4m <- adjacency(
datExpr_inv4m,
power = pow_inv4m,
type = "unsigned"
)
cat("CTRL adjacency dimensions:", dim(adj_ctrl), "\n")
## CTRL adjacency dimensions: 776 776
cat("Inv4m adjacency dimensions:", dim(adj_inv4m), "\n")
## Inv4m adjacency dimensions: 776 776
modules_all_summary <- module_deg_stats_all %>%
group_by(module) %>%
summarise(
n_genes = n(),
n_upregulated = sum(regulation == "Upregulated", na.rm = TRUE),
n_downregulated = sum(regulation == "Downregulated", na.rm = TRUE),
mean_logFC = mean(logFC, na.rm = TRUE),
median_logFC = median(logFC, na.rm = TRUE),
n_in_inv4m = sum(in_Inv4m, na.rm = TRUE),
n_in_cis = sum(in_cis, na.rm = TRUE),
pct_in_cis = 100 * n_in_cis / n_genes,
.groups = "drop"
) %>%
arrange(desc(abs(mean_logFC)))
cat("Module summary (all samples):\n")
## Module summary (all samples):
print(modules_all_summary)
## # A tibble: 14 × 9
## module n_genes n_upregulated n_downregulated mean_logFC median_logFC
## <chr> <int> <int> <int> <dbl> <dbl>
## 1 ME14 10 1 9 -0.711 -0.370
## 2 ME8 30 28 2 0.588 0.548
## 3 ME9 29 8 21 -0.430 -0.434
## 4 ME7 31 0 31 -0.388 -0.325
## 5 ME10 29 0 29 -0.363 -0.294
## 6 ME3 60 4 56 -0.319 -0.287
## 7 ME2 77 5 72 -0.300 -0.373
## 8 ME4 50 4 46 -0.221 -0.316
## 9 ME1 340 140 200 -0.189 -0.350
## 10 ME12 20 11 9 0.182 0.383
## 11 ME11 21 5 16 -0.168 -0.321
## 12 ME13 13 6 7 -0.155 -0.279
## 13 ME5 35 9 26 -0.155 -0.410
## 14 ME6 31 16 15 -0.00166 0.335
## # ℹ 3 more variables: n_in_inv4m <int>, n_in_cis <int>, pct_in_cis <dbl>
write_csv(modules_all_summary, file.path(paths$intermediate, "modules_all_summary.csv"))
modules_ctrl_summary <- module_deg_stats_ctrl %>%
group_by(module) %>%
summarise(
n_genes = n(),
n_upregulated = sum(regulation == "Upregulated", na.rm = TRUE),
n_downregulated = sum(regulation == "Downregulated", na.rm = TRUE),
mean_logFC = mean(logFC, na.rm = TRUE),
median_logFC = median(logFC, na.rm = TRUE),
n_in_inv4m = sum(in_Inv4m, na.rm = TRUE),
n_in_cis = sum(in_cis, na.rm = TRUE),
pct_in_cis = 100 * n_in_cis / n_genes,
.groups = "drop"
) %>%
arrange(desc(abs(mean_logFC)))
modules_inv4m_summary <- module_deg_stats_inv4m %>%
group_by(module) %>%
summarise(
n_genes = n(),
n_upregulated = sum(regulation == "Upregulated", na.rm = TRUE),
n_downregulated = sum(regulation == "Downregulated", na.rm = TRUE),
mean_logFC = mean(logFC, na.rm = TRUE),
median_logFC = median(logFC, na.rm = TRUE),
n_in_inv4m = sum(in_Inv4m, na.rm = TRUE),
n_in_cis = sum(in_cis, na.rm = TRUE),
pct_in_cis = 100 * n_in_cis / n_genes,
.groups = "drop"
) %>%
arrange(desc(abs(mean_logFC)))
cat("\nModule summary (CTRL):\n")
##
## Module summary (CTRL):
print(modules_ctrl_summary)
## # A tibble: 19 × 9
## module n_genes n_upregulated n_downregulated mean_logFC median_logFC
## <chr> <int> <int> <int> <dbl> <dbl>
## 1 ME16 15 12 3 3.21 4.98
## 2 ME0 27 17 10 0.984 0.763
## 3 ME11 28 7 21 -0.857 -0.399
## 4 ME8 36 12 24 -0.819 -0.428
## 5 ME17 15 4 11 0.752 -0.589
## 6 ME9 36 16 20 -0.702 -0.360
## 7 ME10 30 4 26 -0.683 -0.326
## 8 ME2 102 27 75 -0.655 -0.351
## 9 ME7 39 10 29 -0.524 -0.265
## 10 ME1 116 13 103 -0.509 -0.365
## 11 ME5 56 27 29 0.426 -0.231
## 12 ME18 13 2 11 -0.245 -0.305
## 13 ME4 63 36 27 0.236 0.478
## 14 ME6 40 18 22 -0.178 -0.501
## 15 ME3 82 13 69 -0.173 -0.305
## 16 ME12 27 5 22 -0.115 -0.242
## 17 ME15 16 3 13 -0.0859 -0.376
## 18 ME13 18 7 11 0.0708 -0.251
## 19 ME14 17 4 13 -0.0506 -0.454
## # ℹ 3 more variables: n_in_inv4m <int>, n_in_cis <int>, pct_in_cis <dbl>
cat("\nModule summary (Inv4m):\n")
##
## Module summary (Inv4m):
print(modules_inv4m_summary)
## # A tibble: 19 × 9
## module n_genes n_upregulated n_downregulated mean_logFC median_logFC
## <chr> <int> <int> <int> <dbl> <dbl>
## 1 ME12 16 1 15 -6.32 -6.73
## 2 ME8 26 22 4 1.93 0.730
## 3 ME7 33 18 15 1.50 0.272
## 4 ME18 11 10 1 0.923 0.858
## 5 ME11 17 10 7 0.862 0.534
## 6 ME13 15 0 15 -0.839 -0.312
## 7 ME17 12 10 2 0.689 0.686
## 8 ME6 38 1 37 -0.687 -0.390
## 9 ME16 12 3 9 -0.416 -0.662
## 10 ME10 19 1 18 -0.416 -0.323
## 11 ME3 45 1 44 -0.404 -0.294
## 12 ME0 288 89 199 -0.352 -0.370
## 13 ME2 50 22 28 0.300 -0.199
## 14 ME9 21 3 18 -0.241 -0.374
## 15 ME5 40 13 27 -0.209 -0.361
## 16 ME1 62 14 48 -0.201 -0.315
## 17 ME15 13 4 9 0.155 -0.263
## 18 ME14 14 2 12 -0.0946 -0.350
## 19 ME4 44 13 31 0.0120 -0.390
## # ℹ 3 more variables: n_in_inv4m <int>, n_in_cis <int>, pct_in_cis <dbl>
write_csv(modules_ctrl_summary, file.path(paths$intermediate, "modules_ctrl_summary.csv"))
write_csv(modules_inv4m_summary, file.path(paths$intermediate, "modules_inv4m_summary.csv"))
cat("\nCalculating intramodular connectivity (all samples)...\n")
##
## Calculating intramodular connectivity (all samples)...
connectivity_all <- intramodularConnectivity(adj_all, net_all$colors)
hubs_all <- mod_all %>%
mutate(
module_numeric = as.numeric(gsub("ME", "", module))
) %>%
left_join(
connectivity_all %>%
rownames_to_column("gene") %>%
as_tibble(),
by = "gene"
) %>%
left_join(
effects_df %>%
filter(predictor == "Inv4m") %>%
select(gene, locus_label, desc_merged, logFC, adj.P.Val,
regulation, in_Inv4m, in_cis),
by = "gene"
) %>%
arrange(module, desc(kWithin))
cat("Hub genes calculated (all samples)\n")
## Hub genes calculated (all samples)
cat(" Mean kWithin:", round(mean(hubs_all$kWithin, na.rm = TRUE), 3), "\n")
## Mean kWithin: 15.734
cat(" Median kWithin:", round(median(hubs_all$kWithin, na.rm = TRUE), 3), "\n")
## Median kWithin: 6.15
hubs_all_top <- hubs_all %>%
filter(abs(logFC) >= 1.5) %>%
group_by(module) %>%
arrange(desc(kWithin)) %>%
slice_head(n = 10) %>%
ungroup()
cat("\nTop 10 hub genes per module (all samples):\n")
##
## Top 10 hub genes per module (all samples):
print(hubs_all_top %>%
select(module, gene, locus_label, kWithin, logFC, in_cis))
## # A tibble: 21 × 6
## module gene locus_label kWithin logFC in_cis
## <chr> <chr> <chr> <dbl> <dbl> <lgl>
## 1 ME1 Zm00001eb264460 ychf 89.8 6.54 FALSE
## 2 ME1 Zm00001eb213070 gbp1 89.5 8.97 FALSE
## 3 ME1 Zm00001eb189640 fam32a_2 87.6 5.58 TRUE
## 4 ME1 Zm00001eb192880 tet2 86.8 -6.95 TRUE
## 5 ME1 Zm00001eb192850 <NA> 86.4 -1.99 TRUE
## 6 ME1 Zm00001eb187280 gpm458 85.7 9.91 TRUE
## 7 ME1 Zm00001eb196760 <NA> 85.0 -8.51 TRUE
## 8 ME1 Zm00001eb126050 pip5k1 84.9 4.71 FALSE
## 9 ME1 Zm00001eb190240 <NA> 83.6 -9.33 TRUE
## 10 ME1 Zm00001eb191790 jmj2 82.4 -3.55 TRUE
## # ℹ 11 more rows
write_csv(hubs_all, file.path(paths$intermediate, "hubs_all.csv"))
write_csv(hubs_all_top, file.path(paths$intermediate, "hubs_all_top.csv"))
cat("\nHub gene files exported (all samples)\n")
##
## Hub gene files exported (all samples)
connectivity_threshold_all <- quantile(hubs_all$kWithin, 0.75, na.rm = TRUE)
hubs_all_cis <- hubs_all %>%
filter(in_cis == TRUE, kWithin > connectivity_threshold_all) %>%
arrange(desc(kWithin))
cat("\nCis hub genes (all samples, kWithin > 75th percentile =",
round(connectivity_threshold_all, 2), "):", nrow(hubs_all_cis), "\n")
##
## Cis hub genes (all samples, kWithin > 75th percentile = 16.26 ): 149
if (nrow(hubs_all_cis) > 0) {
print(hubs_all_cis %>%
select(module, gene, locus_label, kWithin, logFC, in_Inv4m))
write_csv(hubs_all_cis, file.path(paths$intermediate, "hubs_all_cis.csv"))
}
## module gene locus_label kWithin logFC in_Inv4m
## 1 ME1 Zm00001eb189640 fam32a_2 87.60824 5.5840769 FALSE
## 2 ME1 Zm00001eb192880 tet2 86.82678 -6.9457261 TRUE
## 3 ME1 Zm00001eb192850 <NA> 86.39492 -1.9946674 TRUE
## 4 ME1 Zm00001eb187280 gpm458 85.67797 9.9097830 FALSE
## 5 ME1 Zm00001eb196760 <NA> 85.01359 -8.5057276 FALSE
## 6 ME1 Zm00001eb190240 <NA> 83.58921 -9.3262051 FALSE
## 7 ME1 Zm00001eb191790 jmj2 82.43619 -3.5479196 TRUE
## 8 ME1 Zm00001eb191820 jmj4 81.50158 -2.8826897 TRUE
## 9 ME1 Zm00001eb192630 <NA> 80.50286 -1.7229795 TRUE
## 10 ME1 Zm00001eb188750 hsp70-17 79.53606 -6.5893888 FALSE
## 11 ME1 Zm00001eb189580 <NA> 78.83744 7.5068431 FALSE
## 12 ME1 Zm00001eb189350 <NA> 78.79253 -9.6708956 FALSE
## 13 ME1 Zm00001eb188420 <NA> 78.02280 -6.8786411 FALSE
## 14 ME1 Zm00001eb188430 <NA> 78.02280 -6.8786411 FALSE
## 15 ME1 Zm00001eb190750 jmj21 77.44843 2.3124566 TRUE
## 16 ME1 Zm00001eb193870 ftshi3 77.25884 -5.4979858 TRUE
## 17 ME1 Zm00001eb193260 trxl2 77.15654 -6.0252641 TRUE
## 18 ME1 Zm00001eb188240 ropgef14 75.78290 -3.0530829 FALSE
## 19 ME1 Zm00001eb190670 zfrvt 75.23342 1.9515558 TRUE
## 20 ME1 Zm00001eb192170 <NA> 75.18992 -7.4688231 TRUE
## 21 ME1 Zm00001eb192970 <NA> 74.48943 -6.3077016 TRUE
## 22 ME1 Zm00001eb192330 roc4-1 73.64992 5.0518680 TRUE
## 23 ME1 Zm00001eb192080 <NA> 73.55259 -1.4839851 TRUE
## 24 ME1 Zm00001eb188070 arf1 73.20311 -2.8277808 FALSE
## 25 ME1 Zm00001eb188460 dof43 72.69711 -5.1335736 FALSE
## 26 ME1 Zm00001eb189200 <NA> 72.62853 -6.1597046 FALSE
## 27 ME1 Zm00001eb189170 <NA> 72.07607 -1.8651869 FALSE
## 28 ME1 Zm00001eb192160 <NA> 72.06273 -6.8776969 TRUE
## 29 ME1 Zm00001eb190350 <NA> 71.92532 -1.8752894 FALSE
## 30 ME1 Zm00001eb193620 zfcchc 71.71367 -8.8405439 TRUE
## 31 ME1 Zm00001eb193960 <NA> 70.54737 -1.2610279 TRUE
## 32 ME1 Zm00001eb193130 <NA> 70.09817 -5.6331437 TRUE
## 33 ME1 Zm00001eb188610 <NA> 69.58071 -5.9771086 FALSE
## 34 ME1 Zm00001eb194200 <NA> 69.01066 -1.3199224 TRUE
## 35 ME1 Zm00001eb188570 his2b5 67.77224 -5.3380982 FALSE
## 36 ME1 Zm00001eb189910 metrs1 67.49021 -2.9121795 FALSE
## 37 ME1 Zm00001eb195370 <NA> 67.29793 1.6383226 FALSE
## 38 ME1 Zm00001eb188330 <NA> 66.21055 -1.1528583 FALSE
## 39 ME1 Zm00001eb189820 <NA> 66.13757 5.2973522 FALSE
## 40 ME1 Zm00001eb195120 <NA> 63.96160 -1.2011818 FALSE
## 41 ME1 Zm00001eb190090 <NA> 63.91537 7.1778814 FALSE
## 42 ME1 Zm00001eb189570 <NA> 63.76418 5.5151847 FALSE
## 43 ME1 Zm00001eb192180 <NA> 63.52667 -5.1958191 TRUE
## 44 ME1 Zm00001eb196780 <NA> 63.50148 -6.8610167 FALSE
## 45 ME1 Zm00001eb195860 <NA> 62.97926 -1.4718951 FALSE
## 46 ME1 Zm00001eb192980 <NA> 62.90425 -1.1956236 TRUE
## 47 ME1 Zm00001eb187430 his402 62.49404 -8.1694502 FALSE
## 48 ME1 Zm00001eb196600 r3h1 62.36288 -2.2746235 FALSE
## 49 ME1 Zm00001eb195800 <NA> 61.61740 1.6858934 FALSE
## 50 ME1 Zm00001eb196830 <NA> 58.89420 -5.3301933 FALSE
## 51 ME1 Zm00001eb188220 <NA> 57.93621 -2.9845203 FALSE
## 52 ME1 Zm00001eb191990 pig3 56.96594 1.9529439 TRUE
## 53 ME1 Zm00001eb192900 <NA> 56.82845 -0.8713048 TRUE
## 54 ME1 Zm00001eb187890 smn 56.73194 3.9673651 FALSE
## 55 ME1 Zm00001eb189900 hb115 56.64169 -2.0719433 FALSE
## 56 ME1 Zm00001eb188990 fad1 55.83089 -2.1286584 FALSE
## 57 ME1 Zm00001eb191260 <NA> 55.65987 -1.3868319 TRUE
## 58 ME1 Zm00001eb191320 <NA> 55.17959 -1.2326403 TRUE
## 59 ME1 Zm00001eb194300 paspa 54.89177 -4.0237443 TRUE
## 60 ME1 Zm00001eb194700 <NA> 54.22765 -2.1314797 TRUE
## 61 ME1 Zm00001eb195080 <NA> 54.10796 -4.8397541 FALSE
## 62 ME1 Zm00001eb196200 <NA> 53.91908 -2.3537716 FALSE
## 63 ME1 Zm00001eb196170 mfsd2 53.81298 -1.6469160 FALSE
## 64 ME1 Zm00001eb194380 sec6 53.67194 2.8225461 TRUE
## 65 ME1 Zm00001eb195940 <NA> 53.59613 -0.6912339 FALSE
## 66 ME1 Zm00001eb189390 <NA> 52.86676 -1.4479467 FALSE
## 67 ME1 Zm00001eb195560 <NA> 52.28661 -5.7425424 FALSE
## 68 ME1 Zm00001eb191710 <NA> 52.07647 -0.7872265 TRUE
## 69 ME1 Zm00001eb196880 npi270 51.14548 1.3003241 FALSE
## 70 ME1 Zm00001eb193300 <NA> 50.09204 -4.8842457 TRUE
## 71 ME1 Zm00001eb189700 <NA> 49.47098 4.4551791 FALSE
## 72 ME1 Zm00001eb190500 <NA> 48.73563 -0.5672530 TRUE
## 73 ME1 Zm00001eb189420 <NA> 48.51522 -0.6316545 FALSE
## 74 ME1 Zm00001eb196100 abf2 48.35790 1.7886406 FALSE
## 75 ME1 Zm00001eb196920 <NA> 48.15056 -2.8468789 FALSE
## 76 ME1 Zm00001eb187030 <NA> 47.89935 -3.2839911 FALSE
## 77 ME1 Zm00001eb195790 <NA> 46.91848 3.1218159 FALSE
## 78 ME1 Zm00001eb196240 pin4 46.62636 -3.3660865 FALSE
## 79 ME1 Zm00001eb191430 <NA> 45.97349 -1.0298254 TRUE
## 80 ME1 Zm00001eb189920 aldh2 45.25686 2.6863326 FALSE
## 81 ME1 Zm00001eb192120 <NA> 44.43765 -2.3947086 TRUE
## 82 ME1 Zm00001eb191150 <NA> 44.17096 1.4618638 TRUE
## 83 ME1 Zm00001eb192250 <NA> 44.16916 -0.8130420 TRUE
## 84 ME1 Zm00001eb188030 <NA> 43.46637 4.2180694 FALSE
## 85 ME1 Zm00001eb195970 <NA> 43.25589 -2.6358333 FALSE
## 86 ME1 Zm00001eb189710 <NA> 43.25525 3.7158176 FALSE
## 87 ME1 Zm00001eb193120 <NA> 42.19165 2.3770809 TRUE
## 88 ME1 Zm00001eb187870 <NA> 40.95528 2.5341132 FALSE
## 89 ME1 Zm00001eb195900 ms30 40.14422 -2.8851397 FALSE
## 90 ME1 Zm00001eb196090 pmei47 40.08216 5.3555573 FALSE
## 91 ME1 Zm00001eb191620 <NA> 39.12616 -1.4519423 TRUE
## 92 ME1 Zm00001eb189120 <NA> 39.01994 -0.5983260 FALSE
## 93 ME1 Zm00001eb192580 <NA> 38.89437 2.7379520 TRUE
## 94 ME1 Zm00001eb195600 ptk2 38.69780 -1.6043996 FALSE
## 95 ME1 Zm00001eb190660 <NA> 38.13736 0.8335129 TRUE
## 96 ME1 Zm00001eb192360 <NA> 37.18787 -0.8833142 TRUE
## 97 ME1 Zm00001eb194180 <NA> 35.83976 4.0975260 TRUE
## 98 ME1 Zm00001eb188700 hsp70-7 35.55379 -1.0096395 FALSE
## 99 ME1 Zm00001eb193470 <NA> 34.95628 1.4143723 TRUE
## 100 ME1 Zm00001eb191560 <NA> 34.69698 -1.3020994 TRUE
## 101 ME1 Zm00001eb197020 <NA> 34.14726 -0.9050130 FALSE
## 102 ME1 Zm00001eb192350 <NA> 33.58008 -2.2712141 TRUE
## 103 ME1 Zm00001eb188680 <NA> 33.08431 -0.8451098 FALSE
## 104 ME1 Zm00001eb190880 <NA> 32.27555 0.4458239 TRUE
## 105 ME1 Zm00001eb193560 <NA> 32.26116 -0.4043503 TRUE
## 106 ME1 Zm00001eb186400 <NA> 31.94061 7.7457413 FALSE
## 107 ME1 Zm00001eb196500 <NA> 31.75273 -1.4366628 FALSE
## 108 ME1 Zm00001eb192460 <NA> 31.57491 0.8579075 TRUE
## 109 ME1 Zm00001eb192400 <NA> 30.74269 -0.8670134 TRUE
## 110 ME1 Zm00001eb195050 <NA> 30.11611 -0.3595747 FALSE
## 111 ME1 Zm00001eb195620 <NA> 29.69879 -1.8772725 FALSE
## 112 ME1 Zm00001eb196180 tom2 29.66175 -1.6691084 FALSE
## 113 ME1 Zm00001eb193970 fbl41 29.47838 -1.2902767 TRUE
## 114 ME1 Zm00001eb195290 <NA> 28.37934 1.4354492 FALSE
## 115 ME1 Zm00001eb186390 <NA> 27.82918 5.2882507 FALSE
## 116 ME1 Zm00001eb194240 <NA> 26.86493 0.7327490 TRUE
## 117 ME1 Zm00001eb195830 zhd18 26.85401 -4.3885577 FALSE
## 118 ME1 Zm00001eb186410 <NA> 26.23715 4.9770347 FALSE
## 119 ME1 Zm00001eb196550 <NA> 26.14473 -1.3159081 FALSE
## 120 ME1 Zm00001eb194270 flz22 26.00829 2.6619638 TRUE
## 121 ME1 Zm00001eb190190 <NA> 25.87007 0.5405926 FALSE
## 122 ME1 Zm00001eb195480 <NA> 25.73144 -0.2307286 FALSE
## 123 ME1 Zm00001eb192840 <NA> 25.10500 -1.6572585 TRUE
## 124 ME1 Zm00001eb196630 ufgt4 24.83341 3.9050799 FALSE
## 125 ME1 Zm00001eb189510 dof9 24.39562 0.5734774 FALSE
## 126 ME1 Zm00001eb195180 <NA> 24.29245 -2.3220037 FALSE
## 127 ME1 Zm00001eb188890 <NA> 23.53009 -0.5838285 FALSE
## 128 ME1 Zm00001eb192060 <NA> 23.43976 -0.3905260 TRUE
## 129 ME1 Zm00001eb197300 <NA> 23.13444 -2.6729304 FALSE
## 130 ME1 Zm00001eb193510 <NA> 21.54209 -0.5357527 TRUE
## 131 ME1 Zm00001eb194590 <NA> 20.65390 -0.8451295 TRUE
## 132 ME1 Zm00001eb191680 <NA> 20.39198 -0.3345883 TRUE
## 133 ME1 Zm00001eb188550 rz567b(klc) 20.23608 2.3169023 FALSE
## 134 ME1 Zm00001eb194080 <NA> 20.08124 -4.9094407 TRUE
## 135 ME1 Zm00001eb194610 <NA> 19.37911 -2.7287497 TRUE
## 136 ME1 Zm00001eb195960 <NA> 19.36948 0.7766816 FALSE
## 137 ME1 Zm00001eb192340 vps3 18.61204 -0.6099370 TRUE
## 138 ME1 Zm00001eb194570 <NA> 18.57230 -1.7845376 TRUE
## 139 ME1 Zm00001eb193660 nii2 18.30599 -1.9515935 TRUE
## 140 ME1 Zm00001eb186450 <NA> 18.15185 -0.3732141 FALSE
## 141 ME1 Zm00001eb187460 <NA> 17.96826 2.1857894 FALSE
## 142 ME1 Zm00001eb191780 nactf51 17.67973 -3.5605202 TRUE
## 143 ME1 Zm00001eb188050 <NA> 17.42347 -1.3477292 FALSE
## 144 ME1 Zm00001eb190450 <NA> 17.06172 -1.5667953 FALSE
## 145 ME1 Zm00001eb194600 <NA> 17.02625 -0.8923818 TRUE
## 146 ME1 Zm00001eb191000 <NA> 16.80309 -3.5079674 TRUE
## 147 ME1 Zm00001eb188350 <NA> 16.72521 -0.8678902 FALSE
## 148 ME1 Zm00001eb192260 <NA> 16.55576 -0.4354540 TRUE
## 149 ME1 Zm00001eb194680 rop8 16.38496 -1.6676644 TRUE
cat("\nCalculating intramodular connectivity (genotype-specific)...\n")
##
## Calculating intramodular connectivity (genotype-specific)...
connectivity_ctrl <- intramodularConnectivity(adj_ctrl, net_ctrl$colors)
connectivity_inv4m <- intramodularConnectivity(adj_inv4m, net_inv4m$colors)
hubs_ctrl <- mod_ctrl %>%
mutate(
module_numeric = as.numeric(gsub("ME", "", module))
) %>%
left_join(
connectivity_ctrl %>%
rownames_to_column("gene") %>%
as_tibble(),
by = "gene"
) %>%
left_join(
effects_df %>%
filter(predictor == "Inv4m") %>%
select(gene, locus_label, desc_merged, logFC, adj.P.Val,
regulation, in_Inv4m, in_cis),
by = "gene"
) %>%
arrange(module, desc(kWithin))
hubs_inv4m <- mod_inv4m %>%
mutate(
module_numeric = as.numeric(gsub("ME", "", module))
) %>%
left_join(
connectivity_inv4m %>%
rownames_to_column("gene") %>%
as_tibble(),
by = "gene"
) %>%
left_join(
effects_df %>%
filter(predictor == "Inv4m") %>%
select(gene, locus_label, desc_merged, logFC, adj.P.Val,
regulation, in_Inv4m, in_cis),
by = "gene"
) %>%
arrange(module, desc(kWithin))
cat("Hub genes calculated\n")
## Hub genes calculated
cat(" CTRL - Mean kWithin:", round(mean(hubs_ctrl$kWithin, na.rm = TRUE), 3), "\n")
## CTRL - Mean kWithin: 6.166
cat(" CTRL - Median kWithin:", round(median(hubs_ctrl$kWithin, na.rm = TRUE), 3), "\n")
## CTRL - Median kWithin: 3.742
cat(" Inv4m - Mean kWithin:", round(mean(hubs_inv4m$kWithin, na.rm = TRUE), 3), "\n")
## Inv4m - Mean kWithin: 0.443
cat(" Inv4m - Median kWithin:", round(median(hubs_inv4m$kWithin, na.rm = TRUE), 3), "\n")
## Inv4m - Median kWithin: 0.164
hubs_ctrl_top <- hubs_ctrl %>%
filter(abs(logFC) >= 1.5) %>%
group_by(module) %>%
arrange(desc(kWithin)) %>%
slice_head(n = 10) %>%
ungroup()
hubs_inv4m_top <- hubs_inv4m %>%
filter(abs(logFC) >= 1.5) %>%
group_by(module) %>%
arrange(desc(kWithin)) %>%
slice_head(n = 10) %>%
ungroup()
cat("\nTop 10 hub genes per module (CTRL):\n")
##
## Top 10 hub genes per module (CTRL):
print(hubs_ctrl_top %>%
select(module, gene, locus_label, kWithin, logFC, in_cis))
## # A tibble: 123 × 6
## module gene locus_label kWithin logFC in_cis
## <chr> <chr> <chr> <dbl> <dbl> <lgl>
## 1 ME0 Zm00001eb187440 <NA> 0.174 -2.00 TRUE
## 2 ME0 Zm00001eb166340 <NA> 0.154 -4.62 FALSE
## 3 ME0 Zm00001eb192580 <NA> 0.113 2.74 TRUE
## 4 ME0 Zm00001eb189820 <NA> 0.112 5.30 TRUE
## 5 ME0 Zm00001eb187870 <NA> 0.107 2.53 TRUE
## 6 ME0 Zm00001eb334460 <NA> 0.0673 4.25 FALSE
## 7 ME0 Zm00001eb186840 <NA> 0.0592 1.69 TRUE
## 8 ME0 Zm00001eb189700 <NA> 0.0544 4.46 TRUE
## 9 ME0 Zm00001eb194570 <NA> 0.0419 -1.78 TRUE
## 10 ME0 Zm00001eb257900 mcm5 0.0395 -2.74 FALSE
## # ℹ 113 more rows
cat("\nTop 10 hub genes per module (Inv4m):\n")
##
## Top 10 hub genes per module (Inv4m):
print(hubs_inv4m_top %>%
select(module, gene, locus_label, kWithin, logFC, in_cis))
## # A tibble: 88 × 6
## module gene locus_label kWithin logFC in_cis
## <chr> <chr> <chr> <dbl> <dbl> <lgl>
## 1 ME0 Zm00001eb334460 <NA> 0.0821 4.25 FALSE
## 2 ME0 Zm00001eb371380 <NA> 0.0613 3.36 FALSE
## 3 ME0 Zm00001eb157030 <NA> 0.0577 5.74 FALSE
## 4 ME0 Zm00001eb195370 <NA> 0.0438 1.64 TRUE
## 5 ME0 Zm00001eb197300 <NA> 0.0423 -2.67 TRUE
## 6 ME0 Zm00001eb195620 <NA> 0.0373 -1.88 TRUE
## 7 ME0 Zm00001eb041890 kan1 0.0365 1.61 FALSE
## 8 ME0 Zm00001eb190750 jmj21 0.0296 2.31 TRUE
## 9 ME0 Zm00001eb190450 <NA> 0.0266 -1.57 TRUE
## 10 ME0 Zm00001eb248260 <NA> 0.0264 2.08 FALSE
## # ℹ 78 more rows
write_csv(hubs_ctrl, file.path(paths$intermediate, "hubs_ctrl.csv"))
write_csv(hubs_ctrl_top, file.path(paths$intermediate, "hubs_ctrl_top.csv"))
write_csv(hubs_inv4m, file.path(paths$intermediate, "hubs_inv4m.csv"))
write_csv(hubs_inv4m_top, file.path(paths$intermediate, "hubs_inv4m_top.csv"))
cat("\nHub gene files exported (genotype-specific)\n")
##
## Hub gene files exported (genotype-specific)
connectivity_threshold_ctrl <- quantile(hubs_ctrl$kWithin, 0.75, na.rm = TRUE)
connectivity_threshold_inv4m <- quantile(hubs_inv4m$kWithin, 0.75, na.rm = TRUE)
hubs_ctrl_cis <- hubs_ctrl %>%
filter(in_cis == TRUE, kWithin > connectivity_threshold_ctrl) %>%
arrange(desc(kWithin))
hubs_inv4m_cis <- hubs_inv4m %>%
filter(in_cis == TRUE, kWithin > connectivity_threshold_inv4m) %>%
arrange(desc(kWithin))
cat("\nCis hub genes (CTRL, kWithin > 75th percentile =",
round(connectivity_threshold_ctrl, 2), "):", nrow(hubs_ctrl_cis), "\n")
##
## Cis hub genes (CTRL, kWithin > 75th percentile = 7.68 ): 52
if (nrow(hubs_ctrl_cis) > 0) {
print(hubs_ctrl_cis %>%
select(module, gene, locus_label, kWithin, logFC, in_Inv4m))
write_csv(hubs_ctrl_cis, file.path(paths$intermediate, "hubs_ctrl_cis.csv"))
}
## module gene locus_label kWithin logFC in_Inv4m
## 1 ME1 Zm00001eb192310 <NA> 25.475435 -0.4493898 TRUE
## 2 ME1 Zm00001eb192140 <NA> 22.674525 -0.5957101 TRUE
## 3 ME2 Zm00001eb191790 jmj2 19.354256 -3.5479196 TRUE
## 4 ME3 Zm00001eb195970 <NA> 16.699627 -2.6358333 FALSE
## 5 ME2 Zm00001eb195960 <NA> 15.448279 0.7766816 FALSE
## 6 ME2 Zm00001eb196080 <NA> 15.390397 0.4392978 FALSE
## 7 ME2 Zm00001eb191820 jmj4 15.355796 -2.8826897 TRUE
## 8 ME3 Zm00001eb189030 tola1 15.270774 -0.4167695 FALSE
## 9 ME5 Zm00001eb190360 <NA> 13.943297 -0.9146507 FALSE
## 10 ME5 Zm00001eb194130 rpl29 13.522970 -0.5654438 TRUE
## 11 ME5 Zm00001eb186900 pip1d 13.258877 1.0867126 FALSE
## 12 ME1 Zm00001eb188520 <NA> 13.028288 0.4027907 FALSE
## 13 ME5 Zm00001eb187200 <NA> 12.951922 0.5562427 FALSE
## 14 ME5 Zm00001eb188310 <NA> 12.566939 0.9907331 FALSE
## 15 ME5 Zm00001eb195780 AY104784 12.443780 -0.4335815 FALSE
## 16 ME5 Zm00001eb196030 cipk33 12.380820 -0.7361484 FALSE
## 17 ME2 Zm00001eb191890 ss5 12.354537 -0.6662528 TRUE
## 18 ME2 Zm00001eb195620 <NA> 12.108495 -1.8772725 FALSE
## 19 ME5 Zm00001eb187590 <NA> 11.750940 -1.4486884 FALSE
## 20 ME5 Zm00001eb189210 <NA> 11.533746 -0.3135679 FALSE
## 21 ME6 Zm00001eb193820 <NA> 11.153724 0.4699150 TRUE
## 22 ME2 Zm00001eb194670 <NA> 10.885101 -0.3568029 TRUE
## 23 ME6 Zm00001eb189950 <NA> 10.834141 0.9848448 FALSE
## 24 ME2 Zm00001eb195120 <NA> 10.558953 -1.2011818 FALSE
## 25 ME5 Zm00001eb188220 <NA> 10.199101 -2.9845203 FALSE
## 26 ME5 Zm00001eb189850 <NA> 10.049028 -0.4890672 FALSE
## 27 ME5 Zm00001eb197120 <NA> 9.980317 -0.8027816 FALSE
## 28 ME2 Zm00001eb192360 <NA> 9.853863 -0.8833142 TRUE
## 29 ME6 Zm00001eb192910 lkrsdh1 9.794229 -1.8177127 TRUE
## 30 ME5 Zm00001eb190390 <NA> 9.705346 -0.2950611 FALSE
## 31 ME2 Zm00001eb196310 emp11 9.580118 -0.5777013 FALSE
## 32 ME8 Zm00001eb189140 <NA> 9.337569 -1.1740956 FALSE
## 33 ME8 Zm00001eb190110 <NA> 9.336529 -0.3443905 FALSE
## 34 ME5 Zm00001eb192110 <NA> 9.186318 0.3908126 TRUE
## 35 ME5 Zm00001eb196390 <NA> 9.157441 0.4858115 FALSE
## 36 ME6 Zm00001eb193660 nii2 9.120091 -1.9515935 TRUE
## 37 ME5 Zm00001eb194270 flz22 8.900508 2.6619638 TRUE
## 38 ME4 Zm00001eb193790 tu1 8.812693 -0.4217700 TRUE
## 39 ME6 Zm00001eb189960 <NA> 8.793860 0.9247815 FALSE
## 40 ME5 Zm00001eb187850 <NA> 8.625734 -3.9035202 FALSE
## 41 ME5 Zm00001eb190190 <NA> 8.614837 0.5405926 FALSE
## 42 ME6 Zm00001eb194560 <NA> 8.391949 0.6795678 TRUE
## 43 ME4 Zm00001eb187910 <NA> 8.302450 -0.9566888 FALSE
## 44 ME6 Zm00001eb196880 npi270 8.208947 1.3003241 FALSE
## 45 ME6 Zm00001eb195600 ptk2 8.104428 -1.6043996 FALSE
## 46 ME8 Zm00001eb193040 cpn10 8.018739 0.6148654 TRUE
## 47 ME8 Zm00001eb192960 ereb14 7.988048 -0.4713370 TRUE
## 48 ME5 Zm00001eb189440 <NA> 7.960522 -2.2523055 FALSE
## 49 ME8 Zm00001eb195380 <NA> 7.956684 0.3763682 FALSE
## 50 ME2 Zm00001eb187970 wtf4 7.870608 -0.6071448 FALSE
## 51 ME6 Zm00001eb192640 <NA> 7.784881 0.7591437 TRUE
## 52 ME5 Zm00001eb187240 <NA> 7.741272 -0.2770947 FALSE
cat("\nCis hub genes (Inv4m, kWithin > 75th percentile =",
round(connectivity_threshold_inv4m, 2), "):", nrow(hubs_inv4m_cis), "\n")
##
## Cis hub genes (Inv4m, kWithin > 75th percentile = 0.57 ): 75
if (nrow(hubs_inv4m_cis) > 0) {
print(hubs_inv4m_cis %>%
select(module, gene, locus_label, kWithin, logFC, in_Inv4m))
write_csv(hubs_inv4m_cis, file.path(paths$intermediate, "hubs_inv4m_cis.csv"))
}
## module gene locus_label kWithin logFC in_Inv4m
## 1 ME12 Zm00001eb188610 <NA> 4.5743323 -5.9771086 FALSE
## 2 ME12 Zm00001eb188420 <NA> 4.5743323 -6.8786411 FALSE
## 3 ME12 Zm00001eb188430 <NA> 4.5743323 -6.8786411 FALSE
## 4 ME12 Zm00001eb192880 tet2 4.4522112 -6.9457261 TRUE
## 5 ME1 Zm00001eb192530 <NA> 3.5168297 0.5052828 TRUE
## 6 ME1 Zm00001eb195000 <NA> 3.3884404 -0.3487955 FALSE
## 7 ME7 Zm00001eb186400 <NA> 3.3497706 7.7457413 FALSE
## 8 ME1 Zm00001eb193840 <NA> 3.2291615 -0.5468399 TRUE
## 9 ME7 Zm00001eb186390 <NA> 3.0504198 5.2882507 FALSE
## 10 ME4 Zm00001eb190110 <NA> 2.9558591 -0.3443905 FALSE
## 11 ME1 Zm00001eb191080 <NA> 2.9296961 0.4173646 TRUE
## 12 ME1 Zm00001eb191380 <NA> 2.5665566 -0.5453458 TRUE
## 13 ME1 Zm00001eb191990 pig3 2.5517685 1.9529439 TRUE
## 14 ME12 Zm00001eb188460 dof43 2.5093080 -5.1335736 FALSE
## 15 ME4 Zm00001eb195030 <NA> 2.4144039 4.4482408 FALSE
## 16 ME1 Zm00001eb194950 <NA> 2.4079665 -0.4457779 FALSE
## 17 ME6 Zm00001eb192310 <NA> 2.3318359 -0.4493898 TRUE
## 18 ME6 Zm00001eb188040 <NA> 2.2910721 -0.5558307 FALSE
## 19 ME7 Zm00001eb186770 gras80 2.1872245 1.3572539 FALSE
## 20 ME1 Zm00001eb192900 <NA> 2.0397012 -0.8713048 TRUE
## 21 ME7 Zm00001eb186410 <NA> 1.9831002 4.9770347 FALSE
## 22 ME1 Zm00001eb189120 <NA> 1.9485396 -0.5983260 FALSE
## 23 ME1 Zm00001eb190620 <NA> 1.8098399 0.5535087 TRUE
## 24 ME4 Zm00001eb197320 <NA> 1.7985674 -1.1460387 FALSE
## 25 ME7 Zm00001eb186790 adf7 1.7212075 3.6088351 FALSE
## 26 ME12 Zm00001eb196760 <NA> 1.6342953 -8.5057276 FALSE
## 27 ME3 Zm00001eb189030 tola1 1.5654221 -0.4167695 FALSE
## 28 ME1 Zm00001eb188720 <NA> 1.5510188 -0.2309395 FALSE
## 29 ME4 Zm00001eb193790 tu1 1.5143964 -0.4217700 TRUE
## 30 ME2 Zm00001eb196100 abf2 1.5059538 1.7886406 FALSE
## 31 ME5 Zm00001eb197120 <NA> 1.4792710 -0.8027816 FALSE
## 32 ME1 Zm00001eb188530 <NA> 1.4558796 0.3425578 FALSE
## 33 ME1 Zm00001eb191680 <NA> 1.3680893 -0.3345883 TRUE
## 34 ME1 Zm00001eb192340 vps3 1.3631534 -0.6099370 TRUE
## 35 ME1 Zm00001eb192780 <NA> 1.3350413 -0.4478335 TRUE
## 36 ME4 Zm00001eb190150 <NA> 1.3208462 -1.0198821 FALSE
## 37 ME1 Zm00001eb195280 <NA> 1.3097273 -0.2827958 FALSE
## 38 ME2 Zm00001eb193030 <NA> 1.2897970 0.2651506 TRUE
## 39 ME5 Zm00001eb189140 <NA> 1.2787723 -1.1740956 FALSE
## 40 ME1 Zm00001eb196900 <NA> 1.2683280 1.3223904 FALSE
## 41 ME12 Zm00001eb188570 his2b5 1.2564178 -5.3380982 FALSE
## 42 ME5 Zm00001eb190920 nanmt1 1.1891775 0.3457271 TRUE
## 43 ME1 Zm00001eb196500 <NA> 1.1311407 -1.4366628 FALSE
## 44 ME12 Zm00001eb192170 <NA> 1.1196161 -7.4688231 TRUE
## 45 ME4 Zm00001eb192910 lkrsdh1 1.0435040 -1.8177127 TRUE
## 46 ME2 Zm00001eb193530 <NA> 0.9992925 -0.3054514 TRUE
## 47 ME2 Zm00001eb188310 <NA> 0.9896982 0.9907331 FALSE
## 48 ME7 Zm00001eb186610 <NA> 0.9506719 3.9380702 FALSE
## 49 ME3 Zm00001eb193640 <NA> 0.9489488 0.8931164 TRUE
## 50 ME8 Zm00001eb191370 iaa18 0.9317701 0.3571391 TRUE
## 51 ME7 Zm00001eb189850 <NA> 0.9229845 -0.4890672 FALSE
## 52 ME4 Zm00001eb194330 <NA> 0.9205948 -0.6346064 TRUE
## 53 ME4 Zm00001eb193820 <NA> 0.9203769 0.4699150 TRUE
## 54 ME2 Zm00001eb195340 <NA> 0.9182327 1.0326064 FALSE
## 55 ME4 Zm00001eb190090 <NA> 0.9177928 7.1778814 FALSE
## 56 ME4 Zm00001eb189890 <NA> 0.9084722 -1.0289273 FALSE
## 57 ME1 Zm00001eb196360 <NA> 0.8984926 0.2602678 FALSE
## 58 ME1 Zm00001eb190350 <NA> 0.8981436 -1.8752894 FALSE
## 59 ME11 Zm00001eb187910 <NA> 0.8807480 -0.9566888 FALSE
## 60 ME2 Zm00001eb187200 <NA> 0.8506035 0.5562427 FALSE
## 61 ME2 Zm00001eb196060 <NA> 0.8454635 -0.3284765 FALSE
## 62 ME8 Zm00001eb191020 nactf26 0.8078752 0.2409455 TRUE
## 63 ME6 Zm00001eb192850 <NA> 0.7990997 -1.9946674 TRUE
## 64 ME10 Zm00001eb190360 <NA> 0.7498262 -0.9146507 FALSE
## 65 ME7 Zm00001eb186860 <NA> 0.7472795 3.0258799 FALSE
## 66 ME1 Zm00001eb192060 <NA> 0.7242752 -0.3905260 TRUE
## 67 ME5 Zm00001eb192960 ereb14 0.7241475 -0.4713370 TRUE
## 68 ME11 Zm00001eb196630 ufgt4 0.7214190 3.9050799 FALSE
## 69 ME1 Zm00001eb187890 smn 0.6805194 3.9673651 FALSE
## 70 ME7 Zm00001eb186980 <NA> 0.6516490 0.7492750 FALSE
## 71 ME18 Zm00001eb187460 <NA> 0.6486416 2.1857894 FALSE
## 72 ME4 Zm00001eb195270 bnl10.05 0.6468154 0.6908783 FALSE
## 73 ME7 Zm00001eb189210 <NA> 0.6346461 -0.3135679 FALSE
## 74 ME4 Zm00001eb191500 <NA> 0.5984473 0.3208739 TRUE
## 75 ME11 Zm00001eb196090 pmei47 0.5979830 5.3555573 FALSE
total_cis_all <- sum(module_deg_stats_all$in_cis, na.rm = TRUE)
total_genes_all <- nrow(module_deg_stats_all)
module_cis_enrichment_all <- module_deg_stats_all %>%
group_by(module) %>%
summarise(
n_genes = n(),
n_in_cis = sum(in_cis, na.rm = TRUE),
.groups = "drop"
) %>%
rowwise() %>%
mutate(
fisher_p = {
a <- n_in_cis
b <- n_genes - n_in_cis
c <- total_cis_all - a
d <- (total_genes_all - total_cis_all) - b
contingency <- matrix(c(a, b, c, d), nrow = 2)
fisher.test(contingency)$p.value
},
odds_ratio = {
a <- n_in_cis
b <- n_genes - n_in_cis
c <- total_cis_all - a
d <- (total_genes_all - total_cis_all) - b
(a * d) / (b * c)
}
) %>%
ungroup() %>%
mutate(
fisher_fdr = p.adjust(fisher_p, method = "fdr")
) %>%
arrange(fisher_p)
cat("Cis enrichment analysis (all samples):\n")
## Cis enrichment analysis (all samples):
print(module_cis_enrichment_all)
## # A tibble: 14 × 6
## module n_genes n_in_cis fisher_p odds_ratio fisher_fdr
## <chr> <int> <int> <dbl> <dbl> <dbl>
## 1 ME1 340 213 7.30e-34 6.63 1.02e-32
## 2 ME2 77 3 1.19e-13 0.0546 8.35e-13
## 3 ME3 60 5 5.71e- 8 0.129 2.67e- 7
## 4 ME4 50 4 6.69e- 7 0.126 2.34e- 6
## 5 ME7 31 6 2.40e- 2 0.366 6.73e- 2
## 6 ME10 29 6 5.10e- 2 0.400 1.19e- 1
## 7 ME14 10 1 9.83e- 2 0.173 1.97e- 1
## 8 ME5 35 9 1.13e- 1 0.532 1.98e- 1
## 9 ME9 29 8 2.47e- 1 0.590 3.84e- 1
## 10 ME13 13 3 3.90e- 1 0.468 5.46e- 1
## 11 ME11 21 10 4.97e- 1 1.45 6.32e- 1
## 12 ME12 20 9 6.44e- 1 1.30 7.51e- 1
## 13 ME6 31 13 7.11e- 1 1.15 7.65e- 1
## 14 ME8 30 11 8.51e- 1 0.910 8.51e- 1
cat("\nModules enriched for cis genes (FDR < 0.05):\n")
##
## Modules enriched for cis genes (FDR < 0.05):
cis_enriched_all <- module_cis_enrichment_all %>%
filter(fisher_fdr < 0.05)
if (nrow(cis_enriched_all) > 0) {
print(cis_enriched_all)
} else {
cat(" None\n")
}
## # A tibble: 4 × 6
## module n_genes n_in_cis fisher_p odds_ratio fisher_fdr
## <chr> <int> <int> <dbl> <dbl> <dbl>
## 1 ME1 340 213 7.30e-34 6.63 1.02e-32
## 2 ME2 77 3 1.19e-13 0.0546 8.35e-13
## 3 ME3 60 5 5.71e- 8 0.129 2.67e- 7
## 4 ME4 50 4 6.69e- 7 0.126 2.34e- 6
total_cis_ctrl <- sum(module_deg_stats_ctrl$in_cis, na.rm = TRUE)
total_genes_ctrl <- nrow(module_deg_stats_ctrl)
module_cis_enrichment_ctrl <- module_deg_stats_ctrl %>%
group_by(module) %>%
summarise(
n_genes = n(),
n_in_cis = sum(in_cis, na.rm = TRUE),
.groups = "drop"
) %>%
rowwise() %>%
mutate(
fisher_p = {
a <- n_in_cis
b <- n_genes - n_in_cis
c <- total_cis_ctrl - a
d <- (total_genes_ctrl - total_cis_ctrl) - b
contingency <- matrix(c(a, b, c, d), nrow = 2)
fisher.test(contingency)$p.value
},
odds_ratio = {
a <- n_in_cis
b <- n_genes - n_in_cis
c <- total_cis_ctrl - a
d <- (total_genes_ctrl - total_cis_ctrl) - b
(a * d) / (b * c)
}
) %>%
ungroup() %>%
mutate(
fisher_fdr = p.adjust(fisher_p, method = "fdr")
) %>%
arrange(fisher_p)
total_cis_inv4m <- sum(module_deg_stats_inv4m$in_cis, na.rm = TRUE)
total_genes_inv4m <- nrow(module_deg_stats_inv4m)
module_cis_enrichment_inv4m <- module_deg_stats_inv4m %>%
group_by(module) %>%
summarise(
n_genes = n(),
n_in_cis = sum(in_cis, na.rm = TRUE),
.groups = "drop"
) %>%
rowwise() %>%
mutate(
fisher_p = {
a <- n_in_cis
b <- n_genes - n_in_cis
c <- total_cis_inv4m - a
d <- (total_genes_inv4m - total_cis_inv4m) - b
contingency <- matrix(c(a, b, c, d), nrow = 2)
fisher.test(contingency)$p.value
},
odds_ratio = {
a <- n_in_cis
b <- n_genes - n_in_cis
c <- total_cis_inv4m - a
d <- (total_genes_inv4m - total_cis_inv4m) - b
(a * d) / (b * c)
}
) %>%
ungroup() %>%
mutate(
fisher_fdr = p.adjust(fisher_p, method = "fdr")
) %>%
arrange(fisher_p)
cat("\nCis enrichment analysis (CTRL):\n")
##
## Cis enrichment analysis (CTRL):
print(module_cis_enrichment_ctrl)
## # A tibble: 19 × 6
## module n_genes n_in_cis fisher_p odds_ratio fisher_fdr
## <chr> <int> <int> <dbl> <dbl> <dbl>
## 1 ME1 116 17 0.00000000183 0.227 0.0000000347
## 2 ME5 56 33 0.00166 2.42 0.0158
## 3 ME3 82 20 0.00556 0.474 0.0352
## 4 ME9 36 22 0.00774 2.60 0.0368
## 5 ME6 40 23 0.0186 2.23 0.0707
## 6 ME12 27 5 0.0278 0.348 0.0819
## 7 ME4 63 33 0.0302 1.83 0.0819
## 8 ME0 27 16 0.0423 2.37 0.100
## 9 ME8 36 20 0.0525 2.04 0.111
## 10 ME18 13 2 0.0920 0.282 0.175
## 11 ME11 28 15 0.116 1.86 0.183
## 12 ME15 16 3 0.122 0.358 0.183
## 13 ME10 30 16 0.125 1.85 0.183
## 14 ME14 17 4 0.219 0.479 0.297
## 15 ME7 39 12 0.317 0.689 0.401
## 16 ME16 15 7 0.596 1.39 0.708
## 17 ME2 102 41 0.745 1.07 0.833
## 18 ME17 15 5 0.792 0.785 0.836
## 19 ME13 18 7 1 1.00 1
cat("\nModules enriched for cis genes (CTRL, FDR < 0.05):\n")
##
## Modules enriched for cis genes (CTRL, FDR < 0.05):
cis_enriched_ctrl <- module_cis_enrichment_ctrl %>%
filter(fisher_fdr < 0.05)
if (nrow(cis_enriched_ctrl) > 0) {
print(cis_enriched_ctrl)
} else {
cat(" None\n")
}
## # A tibble: 4 × 6
## module n_genes n_in_cis fisher_p odds_ratio fisher_fdr
## <chr> <int> <int> <dbl> <dbl> <dbl>
## 1 ME1 116 17 0.00000000183 0.227 0.0000000347
## 2 ME5 56 33 0.00166 2.42 0.0158
## 3 ME3 82 20 0.00556 0.474 0.0352
## 4 ME9 36 22 0.00774 2.60 0.0368
cat("\nCis enrichment analysis (Inv4m):\n")
##
## Cis enrichment analysis (Inv4m):
print(module_cis_enrichment_inv4m)
## # A tibble: 19 × 6
## module n_genes n_in_cis fisher_p odds_ratio fisher_fdr
## <chr> <int> <int> <dbl> <dbl> <dbl>
## 1 ME12 16 14 0.0000730 11.5 0.00139
## 2 ME3 45 7 0.000776 0.274 0.00738
## 3 ME1 62 35 0.00405 2.18 0.0257
## 4 ME6 38 7 0.00959 0.341 0.0456
## 5 ME14 14 1 0.0124 0.118 0.0473
## 6 ME4 44 25 0.0161 2.17 0.0508
## 7 ME17 12 1 0.0346 0.141 0.0940
## 8 ME7 33 18 0.0680 1.95 0.161
## 9 ME9 21 5 0.178 0.485 0.346
## 10 ME13 15 3 0.182 0.388 0.346
## 11 ME10 19 5 0.343 0.556 0.564
## 12 ME2 50 16 0.369 0.728 0.564
## 13 ME16 12 3 0.386 0.521 0.564
## 14 ME8 26 12 0.540 1.37 0.732
## 15 ME15 13 6 0.579 1.36 0.734
## 16 ME0 288 115 0.647 1.08 0.769
## 17 ME18 11 5 0.758 1.32 0.847
## 18 ME11 17 7 0.808 1.11 0.852
## 19 ME5 40 16 0.869 1.05 0.869
cat("\nModules enriched for cis genes (Inv4m, FDR < 0.05):\n")
##
## Modules enriched for cis genes (Inv4m, FDR < 0.05):
cis_enriched_inv4m <- module_cis_enrichment_inv4m %>%
filter(fisher_fdr < 0.05)
if (nrow(cis_enriched_inv4m) > 0) {
print(cis_enriched_inv4m)
} else {
cat(" None\n")
}
## # A tibble: 5 × 6
## module n_genes n_in_cis fisher_p odds_ratio fisher_fdr
## <chr> <int> <int> <dbl> <dbl> <dbl>
## 1 ME12 16 14 0.0000730 11.5 0.00139
## 2 ME3 45 7 0.000776 0.274 0.00738
## 3 ME1 62 35 0.00405 2.18 0.0257
## 4 ME6 38 7 0.00959 0.341 0.0456
## 5 ME14 14 1 0.0124 0.118 0.0473
cat("Loading trans co-expression network...\n")
## Loading trans co-expression network...
trans_edges <- read_csv(file.path(paths$intermediate, "inv4m_bipartite_edges.csv"))
cat("Trans network edges loaded:", nrow(trans_edges), "\n")
## Trans network edges loaded: 3766
cat("Columns:", paste(colnames(trans_edges), collapse = ", "), "\n")
## Columns: from, to, r
module_lookup_all <- mod_all %>%
select(gene, module)
trans_edges$in_wgcna_all <- mapply(
function(f, t) {
exists_forward <- any(edge_list_all$from == f & edge_list_all$to == t)
exists_reverse <- any(edge_list_all$from == t & edge_list_all$to == f)
exists_forward | exists_reverse
},
trans_edges$from,
trans_edges$to
)
trans_edges$wgcna_weight_all <- mapply(
function(f, t) {
forward_match <- edge_list_all %>%
filter(from == f & to == t) %>%
pull(weight)
reverse_match <- edge_list_all %>%
filter(from == t & to == f) %>%
pull(weight)
if (length(forward_match) > 0) {
return(forward_match[1])
} else if (length(reverse_match) > 0) {
return(reverse_match[1])
} else {
return(NA_real_)
}
},
trans_edges$from,
trans_edges$to
)
trans_edges_all <- trans_edges %>%
left_join(
module_lookup_all,
by = c("to" = "gene")
) %>%
rename(module_all = module)
cat("\nTrans edges in WGCNA network (all samples):",
sum(trans_edges$in_wgcna_all), "out of", nrow(trans_edges), "\n")
##
## Trans edges in WGCNA network (all samples): 3766 out of 3766
trans_summary_all <- trans_edges_all %>%
filter(!is.na(module_all)) %>%
group_by(module_all) %>%
summarise(
n_trans_edges = n(),
n_in_wgcna = sum(in_wgcna_all, na.rm = TRUE),
pct_in_wgcna = 100 * n_in_wgcna / n_trans_edges,
n_unique_trans_genes = n_distinct(to),
n_unique_cis_regulators = n_distinct(from),
mean_abs_r = mean(abs(r)),
mean_wgcna_weight = mean(wgcna_weight_all, na.rm = TRUE),
.groups = "drop"
) %>%
arrange(desc(n_trans_edges))
cat("\nTrans edges per WGCNA module (all samples):\n")
##
## Trans edges per WGCNA module (all samples):
print(trans_summary_all)
## # A tibble: 5 × 8
## module_all n_trans_edges n_in_wgcna pct_in_wgcna n_unique_trans_genes
## <chr> <int> <int> <dbl> <int>
## 1 ME1 3473 3473 100 33
## 2 ME8 89 89 100 1
## 3 ME14 73 73 100 1
## 4 ME12 68 68 100 1
## 5 ME3 63 63 100 1
## # ℹ 3 more variables: n_unique_cis_regulators <int>, mean_abs_r <dbl>,
## # mean_wgcna_weight <dbl>
write_csv(trans_edges_all, file.path(paths$intermediate, "trans_modules_all.csv"))
write_csv(trans_summary_all, file.path(paths$intermediate, "trans_summary_all.csv"))
cat("\nExported:\n")
##
## Exported:
cat(" trans_modules_all.csv\n")
## trans_modules_all.csv
cat(" trans_summary_all.csv\n")
## trans_summary_all.csv
cat("Calculating TOM (all samples)...\n")
## Calculating TOM (all samples)...
TOM_all <- TOMsimilarity(adj_all)
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.
rownames(TOM_all) <- rownames(adj_all)
colnames(TOM_all) <- colnames(adj_all)
cat("Using WGCNA dendrogram ordering...\n")
## Using WGCNA dendrogram ordering...
all_genes <- colnames(TOM_all)
genes_with_modules <- names(net_all$colors)[net_all$colors != 0]
TOM_all_filtered <- TOM_all[genes_with_modules, genes_with_modules]
if (!is.null(net_all$dendrograms) && length(net_all$dendrograms) > 0) {
full_order <- net_all$dendrograms[[1]]$order
all_genes_ordered <- all_genes[full_order]
genes_in_modules_ordered <- all_genes_ordered[all_genes_ordered %in% genes_with_modules]
gene_order_all <- match(genes_in_modules_ordered, colnames(TOM_all_filtered))
} else {
gene_order_all <- order(net_all$colors[genes_with_modules], genes_with_modules)
}
TOM_all_ordered <- TOM_all_filtered[gene_order_all, gene_order_all]
TOM_all_ordered_rev <- TOM_all_ordered[nrow(TOM_all_ordered):1, ]
gene_names_ordered_all <- colnames(TOM_all_filtered)[gene_order_all]
gene_names_ordered_all_rev <- rev(gene_names_ordered_all)
gene_module_annotation_all <- mod_all %>%
filter(gene %in% gene_names_ordered_all_rev) %>%
select(gene, module) %>%
mutate(gene = factor(gene, levels = gene_names_ordered_all_rev)) %>%
arrange(gene) %>%
column_to_rownames("gene")
pheatmap(
TOM_all_ordered_rev,
annotation_row = gene_module_annotation_all,
annotation_col = gene_module_annotation_all,
cluster_rows = FALSE,
cluster_cols = FALSE,
show_rownames = FALSE,
show_colnames = FALSE,
main = paste0("WGCNA TOM Similarity - All Samples\n(power = ",
pow_all, ")"),
breaks = seq(0, 1, length.out = 101),
color = colorRampPalette(c("white", "red"))(100)
)
write_csv(
as.data.frame(TOM_all_ordered_rev) %>% rownames_to_column("gene"),
file.path(paths$intermediate, "tom_all.csv")
)
cat("\nWGCNA TOM matrix exported (all samples)\n")
##
## WGCNA TOM matrix exported (all samples)
cat("Calculating TOM (genotype-specific)...\n")
## Calculating TOM (genotype-specific)...
TOM_ctrl <- TOMsimilarity(adj_ctrl)
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.
rownames(TOM_ctrl) <- rownames(adj_ctrl)
colnames(TOM_ctrl) <- colnames(adj_ctrl)
TOM_inv4m <- TOMsimilarity(adj_inv4m)
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.
rownames(TOM_inv4m) <- rownames(adj_inv4m)
colnames(TOM_inv4m) <- colnames(adj_inv4m)
cat("TOM matrices calculated\n")
## TOM matrices calculated
all_genes_inv4m <- colnames(TOM_inv4m)
genes_with_modules_inv4m <- names(net_inv4m$colors)[net_inv4m$colors != 0]
TOM_inv4m_filtered <- TOM_inv4m[genes_with_modules_inv4m, genes_with_modules_inv4m]
TOM_ctrl_filtered <- TOM_ctrl[genes_with_modules_inv4m, genes_with_modules_inv4m]
cat("Using Inv4m WGCNA dendrogram ordering...\n")
## Using Inv4m WGCNA dendrogram ordering...
if (!is.null(net_inv4m$dendrograms) && length(net_inv4m$dendrograms) > 0) {
full_order_inv4m <- net_inv4m$dendrograms[[1]]$order
all_genes_ordered_inv4m <- all_genes_inv4m[full_order_inv4m]
genes_in_modules_ordered_inv4m <- all_genes_ordered_inv4m[all_genes_ordered_inv4m %in% genes_with_modules_inv4m]
gene_order_inv4m <- match(genes_in_modules_ordered_inv4m, colnames(TOM_inv4m_filtered))
} else {
gene_order_inv4m <- order(net_inv4m$colors[genes_with_modules_inv4m], genes_with_modules_inv4m)
}
TOM_ctrl_ordered_inv4m <- TOM_ctrl_filtered[gene_order_inv4m, gene_order_inv4m]
TOM_inv4m_ordered_inv4m <- TOM_inv4m_filtered[gene_order_inv4m, gene_order_inv4m]
combined_tom_inv4m <- TOM_ctrl_ordered_inv4m
combined_tom_inv4m[lower.tri(combined_tom_inv4m)] <-
TOM_inv4m_ordered_inv4m[lower.tri(TOM_inv4m_ordered_inv4m)]
combined_tom_inv4m_t <- t(combined_tom_inv4m)
combined_tom_inv4m_rev <- combined_tom_inv4m_t[nrow(combined_tom_inv4m_t):1, ]
gene_names_ordered_inv4m <- colnames(TOM_inv4m_filtered)[gene_order_inv4m]
gene_names_ordered_inv4m_rev <- rev(gene_names_ordered_inv4m)
gene_module_annotation_inv4m <- mod_inv4m %>%
filter(gene %in% gene_names_ordered_inv4m_rev) %>%
select(gene, module) %>%
mutate(gene = factor(gene, levels = gene_names_ordered_inv4m_rev)) %>%
arrange(gene) %>%
column_to_rownames("gene")
pheatmap(
combined_tom_inv4m_rev,
annotation_row = gene_module_annotation_inv4m,
annotation_col = gene_module_annotation_inv4m,
cluster_rows = FALSE,
cluster_cols = FALSE,
show_rownames = FALSE,
show_colnames = FALSE,
main = paste0("WGCNA TOM Similarity (Upper: CTRL [pow=", pow_ctrl,
"], Lower: Inv4m [pow=", pow_inv4m, "])\nClustered by Inv4m"),
breaks = seq(0, 1, length.out = 101),
color = colorRampPalette(c("white", "red"))(100)
)
write_csv(
as.data.frame(TOM_ctrl_ordered_inv4m) %>% rownames_to_column("gene"),
file.path(paths$intermediate, "tom_ctrl_inv4m_order.csv")
)
write_csv(
as.data.frame(TOM_inv4m_ordered_inv4m) %>% rownames_to_column("gene"),
file.path(paths$intermediate, "tom_inv4m_inv4m_order.csv")
)
all_genes_ctrl <- colnames(TOM_ctrl)
genes_with_modules_ctrl <- names(net_ctrl$colors)[net_ctrl$colors != 0]
TOM_ctrl_filtered_ctrl <- TOM_ctrl[genes_with_modules_ctrl, genes_with_modules_ctrl]
TOM_inv4m_filtered_ctrl <- TOM_inv4m[genes_with_modules_ctrl, genes_with_modules_ctrl]
cat("Using CTRL WGCNA dendrogram ordering...\n")
## Using CTRL WGCNA dendrogram ordering...
if (!is.null(net_ctrl$dendrograms) && length(net_ctrl$dendrograms) > 0) {
full_order_ctrl <- net_ctrl$dendrograms[[1]]$order
all_genes_ordered_ctrl <- all_genes_ctrl[full_order_ctrl]
genes_in_modules_ordered_ctrl <- all_genes_ordered_ctrl[all_genes_ordered_ctrl %in% genes_with_modules_ctrl]
gene_order_ctrl <- match(genes_in_modules_ordered_ctrl, colnames(TOM_ctrl_filtered_ctrl))
} else {
gene_order_ctrl <- order(net_ctrl$colors[genes_with_modules_ctrl], genes_with_modules_ctrl)
}
TOM_ctrl_ordered_ctrl <- TOM_ctrl_filtered_ctrl[gene_order_ctrl, gene_order_ctrl]
TOM_inv4m_ordered_ctrl <- TOM_inv4m_filtered_ctrl[gene_order_ctrl, gene_order_ctrl]
combined_tom_ctrl <- TOM_inv4m_ordered_ctrl
combined_tom_ctrl[lower.tri(combined_tom_ctrl)] <-
TOM_ctrl_ordered_ctrl[lower.tri(TOM_ctrl_ordered_ctrl)]
combined_tom_ctrl_t <- t(combined_tom_ctrl)
combined_tom_ctrl_rev <- combined_tom_ctrl_t[nrow(combined_tom_ctrl_t):1, ]
gene_names_ordered_ctrl <- colnames(TOM_ctrl_filtered_ctrl)[gene_order_ctrl]
gene_names_ordered_ctrl_rev <- rev(gene_names_ordered_ctrl)
gene_module_annotation_ctrl <- mod_ctrl %>%
filter(gene %in% gene_names_ordered_ctrl_rev) %>%
select(gene, module) %>%
mutate(gene = factor(gene, levels = gene_names_ordered_ctrl_rev)) %>%
arrange(gene) %>%
column_to_rownames("gene")
pheatmap(
combined_tom_ctrl_rev,
annotation_row = gene_module_annotation_ctrl,
annotation_col = gene_module_annotation_ctrl,
cluster_rows = FALSE,
cluster_cols = FALSE,
show_rownames = FALSE,
show_colnames = FALSE,
main = paste0("WGCNA TOM Similarity (Upper: Inv4m [pow=", pow_inv4m,
"], Lower: CTRL [pow=", pow_ctrl, "])\nClustered by CTRL"),
breaks = seq(0, 1, length.out = 101),
color = colorRampPalette(c("white", "red"))(100)
)
write_csv(
as.data.frame(TOM_ctrl_ordered_ctrl) %>% rownames_to_column("gene"),
file.path(paths$intermediate, "tom_ctrl_ctrl_order.csv")
)
write_csv(
as.data.frame(TOM_inv4m_ordered_ctrl) %>% rownames_to_column("gene"),
file.path(paths$intermediate, "tom_inv4m_ctrl_order.csv")
)
cat("\nGenotype-specific TOM matrices exported\n")
##
## Genotype-specific TOM matrices exported